function [eigen_value] = myEigen(M)

m = trace(M)/3;
K = M-m*eye(3);
q = det(K)/2;
 
p = 0;
for i=1:3
    for j=1:3
        p = p + K(i,j)^2;
    end
end
p = p/6;
 
 
phi = 1/3*acos(q/p^(3/2));
 
% NOTE: the follow formula assume accurate computation and therefor q/p^(3/2) should be in range of [1,-1], 
% but in real code, because of numerical errors, it must be checked. Thus, in case abs(q) >= abs(p^(3/2)), set phi = 0;
if(abs(q) >= abs(p^(3/2)))
    phi = 0;
end
 
if(phi<0)
    phi=phi+pi/3;
end

eigen(1) = m + 2*sqrt(p)*cos(phi);
eigen(2) = m - sqrt(p)*(cos(phi) + sqrt(3)*sin(phi));
eigen(3) = m - sqrt(p)*(cos(phi) - sqrt(3)*sin(phi));

sort(eigen);

eigen_value = zeros(3);
eigen_value(1,1) = eigen(1);
eigen_value(2,2) = eigen(2);
eigen_value(3,3) = eigen(3);

end